PatchSizePilot

Biomass

Aim

Here I study how biomass density changes across treatments in the PatchSizePilot. In particular, I’m studying how biomass density changes across:

ecosystems of different size (this is the case in nature, but we need to make sure that it’s the case also in our experiment)
ecosystems that are connected (through the flow of nutrients) to an ecosystem of the same size or to an ecosystem that is larger
meta-ecosystems in which patch size is the same or in which one patch has most of the area (we keep the total area of the meta-ecosystem constant)

Data manipulation

Import

culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "population", "t0.RData")); t0 = pop_output
load(here("data", "population", "t1.RData")); t1 = pop_output
load(here("data", "population", "t2.RData")); t2 = pop_output
load(here("data", "population", "t3.RData")); t3 = pop_output
load(here("data", "population", "t4.RData")); t4 = pop_output
load(here("data", "population", "t5.RData")); t5 = pop_output
load(here("data", "population", "t6.RData")); t6 = pop_output
load(here("data", "population", "t7.RData")); t7 = pop_output
rm(pop_output)

Tidy

#Column: time
t0$time = NA
t1$time = NA

#Column: replicate_video
t0$replicate_video = 1:12 #In t1 I took 12 videos of a single 
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>%
  rename(replicate_video = replicate)
t7 = t7 %>%
  rename(replicate_video = replicate)

#Create an elongated version of t0 so that each of the 110 cultures can have 12 video replicates at t0.
elongating_t0 = NULL
for (video in 1:nrow(t0)){
  
  for (ID in 1:nrow(culture_info)) {
    
    elongating_t0 = rbind(elongating_t0, t0[video,])
    
    }

  }

ID_vector = rep(1:nrow(culture_info), 
                times = nrow(t0))

elongating_t0$culture_ID = ID_vector

#Merge previous data-sets
t0 = merge(culture_info,elongating_t0, by="culture_ID")
t1 = merge(culture_info,t1, by = "culture_ID")
t2 = merge(culture_info,t2, by = "culture_ID")
t3 = merge(culture_info,t3, by = "culture_ID")
t4 = merge(culture_info,t4, by = "culture_ID")
t5 = merge(culture_info,t5, by = "culture_ID")
t6 = merge(culture_info,t6, by = "culture_ID")
t7 = merge(culture_info,t7, by = "culture_ID")
ds = rbind(t0, t1, t2, t3, t4, t5, t6, t7)
rm(elongating_t0, t0, t1, t2, t3, t4, t5, t6, t7)

ds$time_point[ds$time_point=="t0"] = 0
ds$time_point[ds$time_point=="t1"] = 1
ds$time_point[ds$time_point=="t2"] = 2
ds$time_point[ds$time_point=="t3"] = 3
ds$time_point[ds$time_point=="t4"] = 4
ds$time_point[ds$time_point=="t5"] = 5
ds$time_point[ds$time_point=="t6"] = 6
ds$time_point[ds$time_point=="t7"] = 7
ds$time_point = as.character(ds$time_point)

ds$day = NA
ds$day[ds$time_point== 0] = 0
ds$day[ds$time_point== 1] = 4
ds$day[ds$time_point== 2] = 8
ds$day[ds$time_point== 3] = 12
ds$day[ds$time_point== 4] = 16
ds$day[ds$time_point== 5] = 20
ds$day[ds$time_point== 6] = 24
ds$day[ds$time_point== 7] = 28

#Column: eco_metaeco_type
ds$eco_metaeco_type = factor(ds$eco_metaeco_type, 
                             levels = c('S', 
                                        'S (S_S)', 
                                        'S (S_L)', 
                                        'M', 
                                        'M (M_M)', 
                                        'L', 
                                        'L (L_L)', 
                                        'L (S_L)'))

ecosystems_to_take_off = 60 #Culture number 60 because it was spilled
ds = ds %>%
  filter(! culture_ID %in% ecosystems_to_take_off)

ds_for_evaporation = ds

ds = ds %>% 
  select(culture_ID, 
         patch_size, 
         disturbance, 
         metaecosystem_type, 
         bioarea_per_volume, 
         replicate_video, 
         time_point,
         day,
         metaecosystem, 
         system_nr, 
         eco_metaeco_type)

ds = ds[, c("culture_ID", 
            "system_nr", 
            "disturbance", 
            "time_point",
            "day", 
            "patch_size", 
            "metaecosystem", 
            "metaecosystem_type", 
            "eco_metaeco_type", 
            "replicate_video",
            "bioarea_per_volume")]

Create regional data set

ds_regional = ds %>%
  filter(metaecosystem == "yes") %>%
  group_by(culture_ID, 
           system_nr, 
           disturbance, 
           day, 
           time_point, 
           patch_size, 
           metaecosystem_type) %>%
  summarise(patch_mean_bioarea_across_videos = mean(bioarea_per_volume)) %>%
  group_by(system_nr, disturbance, day, time_point, metaecosystem_type) %>%
  summarise(regional_mean_bioarea = mean(patch_mean_bioarea_across_videos))

metaecosystems_to_take_off = 40 #System 40 was the system of culture 60 that I spilled
ds_regional = ds_regional %>%
  filter(! system_nr %in% metaecosystems_to_take_off)

Medium-Medium vs Small-Large

We want to understand if the regional biomass produced by an ecosystem with a small and a large patch (metaecosystem_type = S_L) is lower than the regional biomass produced by an ecosystem with two medium patches (metaecosystem_type = M_M).

Plot

Let’s start by plotting the single ecosystems to see that everything is fine. To make the patterns clear let’s plot the low disturbance and high disturbance in two different figures.

Let’s then show these data as boxplots.

We can see that the regional biomass was higher for meta-ecosystems with the same size, regardless of disturbance. As both disturbance levels showed the same pattern, I would just keep the one with low disturbance in the publication, as it shows a stronger pattern.

Time as random effect

Tidy

First of all, let’s modify the data set including the regional biomass of our meta-ecosystems. In this data set, we want to have the regional biomass of the meta-ecosystems (averaged first across videos and then across patches) in which we:

Include only the meta-ecosystems in which patches had both medium size (metaecosystem_type = M_M) and meta-ecosystems in which patches had one a small size and the other large size (metaecosystem_type = S_L).
Take off the first two point (day 0 and day = 4). This is because the first perturbation happened only at day 5.

ds_regional_MM_SL_t2t7 = ds_regional %>%
    filter (metaecosystem_type == "M_M" | metaecosystem_type == "S_L", 
            time_point >= 2)

Model selection

Let’s start from the largest mixed effect model makes sense to construct. This will include:

As fixed effect: metaecosystem type (M), disturbance (D), and their interaction (M * D)
As random effect (random intercept and slope): day and the number of the metaecosystem (system_nr).

We are actually going to construct two full models. One with the correlated and one with the uncorrelated random slopes and intercept, the other without

(1) Should we keep the correlation between intercept and slopes (|)?

full_model_correlated = lmer(regional_mean_bioarea ~ 
                             metaecosystem_type  + 
                             disturbance + 
                             metaecosystem_type * disturbance + 
                             (metaecosystem_type | day) + 
                             (disturbance | day) + 
                             (metaecosystem_type*disturbance  | day) +
                             (1 | system_nr) , 
                             data = ds_regional_MM_SL_t2t7, 
                             REML = FALSE)

full_model_uncorrelated = lmer(regional_mean_bioarea ~ 
                    metaecosystem_type  + 
                    disturbance + 
                    metaecosystem_type * disturbance + 
                    (metaecosystem_type || day) + 
                    (disturbance || day) + 
                    (metaecosystem_type*disturbance  || day) +
                    (1 | system_nr) ,
                  data = ds_regional_MM_SL_t2t7, 
                  REML = FALSE)

anova(full_model_uncorrelated, full_model_correlated)

## Data: ds_regional_MM_SL_t2t7
## Models:
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
## full_model_uncorrelated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + ((1 | day) + (0 + metaecosystem_type | day)) + ((1 | day) + (0 + disturbance | day)) + ((1 | day) + (0 + metaecosystem_type | day) + (0 + disturbance | day) + (0 + metaecosystem_type:disturbance | day)) + (1 | system_nr)
##                         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## full_model_correlated     22 1881.2 1942.5 -918.58   1837.2                    
## full_model_uncorrelated   31 1899.3 1985.7 -918.63   1837.3     0  9          1

Yes (according to BIC and AIC). Interestingly enough, the p-value is 1 (how should I interpret this?).

(2) Should we keep M * D ?

no_interaction_model = lmer(regional_mean_bioarea ~ 
                              metaecosystem_type + 
                              disturbance  + 
                              (metaecosystem_type | day) + 
                              (disturbance | day) +
                              (1 | system_nr), 
                            data = ds_regional_MM_SL_t2t7, 
                            REML = FALSE)

anova(full_model_correlated, no_interaction_model)

## Data: ds_regional_MM_SL_t2t7
## Models:
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
##                       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## no_interaction_model    11 1861.7 1892.4 -919.86   1839.7                     
## full_model_correlated   22 1881.2 1942.5 -918.58   1837.2 2.5592 11     0.9954

No (according to AIC and BIC). The p-value is also here almost 1 (how should I interpret this?).

(3) Should we keep (M | day)?

no_metaeco_slopes_model = lmer(regional_mean_bioarea ~ 
                                 metaecosystem_type + 
                                 disturbance  + 
                                 (disturbance | day) +
                                 (1 | system_nr), 
                               data = ds_regional_MM_SL_t2t7, 
                               REML = FALSE)

anova(no_interaction_model, no_metaeco_slopes_model)

## Data: ds_regional_MM_SL_t2t7
## Models:
## no_metaeco_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## no_metaeco_slopes_model    8 1871.0 1893.2 -927.48   1855.0          
## no_interaction_model      11 1861.7 1892.4 -919.86   1839.7 15.239  3
##                         Pr(>Chisq)   
## no_metaeco_slopes_model              
## no_interaction_model      0.001623 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yes (according to AIC and BIC - however, BIC is only slightly lower). The p-value is really low (0.002).

(4) Should we keep (D | day)?

no_disturbance_slopes_model = lmer(regional_mean_bioarea ~ 
                                     metaecosystem_type + 
                                     disturbance  + 
                                     (metaecosystem_type | day) +
                                     (1 | system_nr), 
                                   data = ds_regional_MM_SL_t2t7, 
                                   REML = FALSE)

anova(no_interaction_model, no_disturbance_slopes_model)

## Data: ds_regional_MM_SL_t2t7
## Models:
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
##                             npar    AIC    BIC  logLik deviance  Chisq Df
## no_disturbance_slopes_model    8 1860.6 1882.9 -922.31   1844.6          
## no_interaction_model          11 1861.7 1892.4 -919.86   1839.7 4.9107  3
##                             Pr(>Chisq)
## no_disturbance_slopes_model           
## no_interaction_model            0.1785

No (according to AIC and BIC - however, AIC is only slightly lower). The p-value is low (0.1785).

(5) Should we keep (1 | system_nr)?

no_system_nr = lmer(regional_mean_bioarea ~ 
                      metaecosystem_type + 
                      disturbance  + 
                      (metaecosystem_type | day), 
                    data = ds_regional_MM_SL_t2t7, 
                    REML = FALSE)

anova(no_disturbance_slopes_model, no_system_nr)

## Data: ds_regional_MM_SL_t2t7
## Models:
## no_system_nr: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day)
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
##                             npar    AIC    BIC  logLik deviance Chisq Df
## no_system_nr                   7 1861.5 1881.0 -923.75   1847.5         
## no_disturbance_slopes_model    8 1860.6 1882.9 -922.31   1844.6 2.881  1
##                             Pr(>Chisq)  
## no_system_nr                            
## no_disturbance_slopes_model    0.08963 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yes according to AIC (but no according to BIC). The p-value is low (0.09).

Should we keep these parameters?

Parameter	AIC	BIC
\| (correlation)	Yes	Yes
M * D	No	No
(M \| t)	Yes	Yes (but slightly lower)
(D \| t)	No (but slightly lower)	No
(1 \| ID)	Yes	No

According to AIC, the best model should be M + D + (1 | system_nr). BIC would take out (1 | system_nr). But as we are going for AIC here, our best model should include it.

To make sure everything is okay, let’s take a look at all our models.

anova(full_model_correlated, full_model_uncorrelated, no_interaction_model, no_metaeco_slopes_model, no_disturbance_slopes_model, no_system_nr)

## Data: ds_regional_MM_SL_t2t7
## Models:
## no_system_nr: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day)
## no_metaeco_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_disturbance_slopes_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (1 | system_nr)
## no_interaction_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## full_model_correlated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type * disturbance | day) + (1 | system_nr)
## full_model_uncorrelated: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type * disturbance + ((1 | day) + (0 + metaecosystem_type | day)) + ((1 | day) + (0 + disturbance | day)) + ((1 | day) + (0 + metaecosystem_type | day) + (0 + disturbance | day) + (0 + metaecosystem_type:disturbance | day)) + (1 | system_nr)
##                             npar    AIC    BIC  logLik deviance   Chisq Df
## no_system_nr                   7 1861.5 1881.0 -923.75   1847.5           
## no_metaeco_slopes_model        8 1871.0 1893.2 -927.48   1855.0  0.0000  1
## no_disturbance_slopes_model    8 1860.6 1882.9 -922.31   1844.6 10.3288  0
## no_interaction_model          11 1861.7 1892.4 -919.86   1839.7  4.9107  3
## full_model_correlated         22 1881.2 1942.5 -918.58   1837.2  2.5592 11
## full_model_uncorrelated       31 1899.3 1985.7 -918.63   1837.3  0.0000  9
##                             Pr(>Chisq)
## no_system_nr                          
## no_metaeco_slopes_model         1.0000
## no_disturbance_slopes_model           
## no_interaction_model            0.1785
## full_model_correlated           0.9954
## full_model_uncorrelated         1.0000

Odd, in some of the models chi squared is zero and in some the degrees of freedom are zero.

The best model according to AIC should be M + D + (M | day) + (1 | system number) (no_disturbance_slopes_model). However, it has zero degrees of freedom.

So, the models I’m going to consider as possible models are the following ones:

M + D + (M | t) + (1 | ID) [Best model according to AIC]
M + D + (M | t) + (D | t) + (1 | ID) [Best model according to AIC with (D | t), which had a slightly higher AIC]
M + D + (M | t) [Best model according to BIC]
M + D [Best model according to BIC without (M | t), which had only a slightly higher BIC]

Time fitting

Here we want to fit how biomass changes across time to a function. The biomass of our meta-ecosystems looks like this.

To fit these data, we need to produce (and parameterise) a function that resemble how the biomass first increases and then decreases.

Hank produced the following function:

\[biomass = a_4 * (day - a_5) * e^{a_1(day - a_5)}\]

If we parameterise the function and then plot, it looks as follows.

a1 = -0.1
a4 = 1200
a5 = -1

day = seq(0, 30, 0.01)
biomass = a4*(day-a5) * exp(a1*(day-a5))
plot(biomass ~ day)

Now, let’s find the best parameters (a1, a4, a5) that fit our data.

model$m$getPars()

##           a1           a4           a5 
##   -0.1336182 1216.4268967   -1.6147345

And now let’s plot the function to see how it fits our data.

Let’s then include our predictions into the data set.

ds_regional_predicted_shrunk_type = ds_regional %>%
  mutate(predicted_from_time = a4*(day-a5)*exp(a1*(day-a5))) %>%
  filter(metaecosystem_type == "S_L" | metaecosystem_type == "M_M")

Time as fixed effect

Let’s now work using the fitted data that we found in the previous section.

Tidy

ds_regional_predicted_shrunk_type_n_day = ds_regional_predicted_shrunk_type %>%
  filter(time_point >= 2)

Model selection

Let’s start from the full model.

model_slopes_correlated = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            predicted_from_time * metaecosystem_type +
                            predicted_from_time * disturbance +
                            metaecosystem_type * disturbance +
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time | system_nr),
                        data = ds_regional_predicted_shrunk_type_n_day,
                        REML = FALSE)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

(1) Should we keep the correlation in (T | ID)?

model_slopes_uncorrelated = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            predicted_from_time * metaecosystem_type +
                            predicted_from_time * disturbance +
                            metaecosystem_type * disturbance +
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time || system_nr),
                        data = ds_regional_predicted_shrunk_type_n_day,
                        REML = FALSE,
                        control = lmerControl(optimizer = "Nelder_Mead"))

## Warning: Some predictor variables are on very different scales: consider
## rescaling

anova(model_slopes_uncorrelated, model_slopes_correlated)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_slopes_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
##                           npar    AIC    BIC  logLik deviance  Chisq Df
## model_slopes_uncorrelated   11 1841.8 1872.5 -909.90   1819.8          
## model_slopes_correlated     12 1842.7 1876.2 -909.35   1818.7 1.1016  1
##                           Pr(>Chisq)
## model_slopes_uncorrelated           
## model_slopes_correlated       0.2939

No (according to AIC and BIC - however, AIC is only slightly lower).

(2) Should we keep T * M?

model_no_TM = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            predicted_from_time * disturbance +
                            metaecosystem_type * disturbance +
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time || system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead")
               )

## Warning: Some predictor variables are on very different scales: consider
## rescaling

anova(model_slopes_uncorrelated, model_no_TM)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_TM: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
##                           npar    AIC    BIC logLik deviance Chisq Df
## model_slopes_uncorrelated   11 1841.8 1872.5 -909.9   1819.8         
## model_no_TM                 11 1841.8 1872.5 -909.9   1819.8     0  0
##                           Pr(>Chisq)
## model_slopes_uncorrelated           
## model_no_TM

This is really wired. Df = 0. Let’s then say no then.

(3) Should we keep T * D?

model_no_TD = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            metaecosystem_type * disturbance +
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time || system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead"))

## Warning: Some predictor variables are on very different scales: consider
## rescaling

anova(model_no_TM, model_no_TD)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TM: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_TD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
##             npar    AIC    BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_TM   11 1841.8 1872.5 -909.9   1819.8                    
## model_no_TD   11 1841.8 1872.5 -909.9   1819.8     0  0

This is really wired. Df = 0. Let’s then say no then.

(4) Should we keep M * D?

model_no_MD = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time || system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead"))

## Warning: Some predictor variables are on very different scales: consider
## rescaling

anova(model_no_TD, model_no_MD)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
##             npar    AIC    BIC logLik deviance Chisq Df Pr(>Chisq)
## model_no_TD   11 1841.8 1872.5 -909.9   1819.8                    
## model_no_MD   11 1841.8 1872.5 -909.9   1819.8     0  0

This is really wired. Df = 0. Let’s then say no then.

(5) Should we keep T * M * D?

model_no_TMD = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            (predicted_from_time || system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead"))

anova(model_no_MD, model_no_TMD)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_TMD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
##              npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## model_no_TMD    7 1844.3 1863.8 -915.15   1830.3                      
## model_no_MD    11 1841.8 1872.5 -909.90   1819.8 10.49  4    0.03294 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yes according to AIC (but no according to BIC).

(6) Should we keep (T || ID)?

model_no_random_slope = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            (1 | system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead"))

anova(model_no_MD, model_no_random_slope)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_random_slope: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + (1 | system_nr)
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
##                       npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## model_no_random_slope    6 1842.3 1859.0 -915.15   1830.3                      
## model_no_MD             11 1841.8 1872.5 -909.90   1819.8 10.49  5    0.06249 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Yes according to AIC (but no according to BIC).

(7) Should we introduce the correlation (T | ID) ?

model_no_MD_correlated = lmer(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance + 
                            predicted_from_time * metaecosystem_type * disturbance +
                            (predicted_from_time | system_nr),
               data = ds_regional_predicted_shrunk_type_n_day,
               REML = FALSE,
               control = lmerControl (optimizer = "Nelder_Mead"))

## Warning: Some predictor variables are on very different scales: consider
## rescaling

anova(model_no_MD, model_no_MD_correlated)

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## model_no_MD: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_no_MD_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
##                        npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## model_no_MD              11 1841.8 1872.5 -909.90   1819.8                    
## model_no_MD_correlated   12 1841.5 1875.0 -908.76   1817.5 2.286  1     0.1305

Yes according to AIC (but no according to BIC). This is odd as at the beginning the model told us we shouldn’t keep it.

Should we keep these parameters?

Parameter	AIC	BIC
(T \| ID) (correlation)	No (but only slightly lower)	No
T * M	?	?
T * D	?	?
M * D	?	?
T * M * D	Yes	No
(T \|\| ID)	Yes	No

According to the previous table, then BIC should be the highest for a fixed effect model. Let’s see if this is the case.

fixed_model = lm(regional_mean_bioarea ~ 
                            predicted_from_time + 
                            metaecosystem_type  + 
                            disturbance,
               data = ds_regional_predicted_shrunk_type_n_day)

To make sure everything is okay, let’s take a look at all our models.

anova(model_slopes_correlated, model_slopes_uncorrelated, model_no_random_slope, fixed_model) #Only with models df > 0

## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## fixed_model: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance
## model_no_random_slope: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + (1 | system_nr)
## model_slopes_uncorrelated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## model_slopes_correlated: regional_mean_bioarea ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time * metaecosystem_type + predicted_from_time * disturbance + metaecosystem_type * disturbance + predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
##                           npar    AIC    BIC  logLik deviance   Chisq Df
## fixed_model                  5 1841.7 1855.6 -915.84   1831.7           
## model_no_random_slope        6 1842.3 1859.0 -915.15   1830.3  1.3783  1
## model_slopes_uncorrelated   11 1841.8 1872.5 -909.90   1819.8 10.4897  5
## model_slopes_correlated     12 1842.7 1876.2 -909.35   1818.7  1.1016  1
##                           Pr(>Chisq)  
## fixed_model                           
## model_no_random_slope        0.24038  
## model_slopes_uncorrelated    0.06249 .
## model_slopes_correlated      0.29393  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#anova(model_slopes_correlated, model_slopes_uncorrelated, model_no_TM, model_no_TD, model_no_MD, model_no_TMD, model_no_random_slope, fixed_model) #With all the models

The best model according to AIC is T + M + D + (T || ID) (model_no_MD), as we found before. However, all other models seem to be of similar goodness. According to BIC, the best is (fixed_model)

Therefore, the models I will consider are:

T + M + D + T * M * D + (T || ID) [best model according to AIC]
T + M + D [best model according to BIC]
T + M + D + T * M * D + (T | ID) [best model found introducing the correlation at the end]*

*DOESN’t WORK THOUGH (number of observations (=40) <= number of random effects (=40) )

Selected models

Evaporation

We want to know if there was a systematic bias in the evaporation of different treatments (disturbance, patch size) and whether evaporation changed across time. My expectation would be that we would see a difference among the exchanges 2,3 and the exchanges 4,5,6. This is because in exchange 2,3 cultures were microwaved in 15 tubes for 3 minutes and in exchange 4,5,6 cultures were microwaved in 4 tubes for only 1 minute.

Tidy

#Columns: exchange & evaporation
ds_for_evaporation = gather(ds_for_evaporation, 
                            key = exchange, 
                            value = evaporation, 
                            water_add_after_t2:water_add_after_t6)
ds_for_evaporation[ds_for_evaporation == "water_add_after_t2"] = "2"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t3"] = "3"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t4"] = "4"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t5"] = "5"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t6"] = "6"
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] / 2 #This is because exchange contained the topping up of two exchanges
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] + 2 #We need to add 2 ml to the evaporation that happened at the exchange events 1 and 2. This is because we already added 1 ml of water at exchange 1 and 1 ml of water at exchange 2. 

#Column: nr_of_tubes_in_rack
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 1] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 2] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 3] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 4] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 5] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 6] = 4

Plot

It seems like there is no real difference across time, disturbance, or patch type. However, we could also run a mixed effect model to show that they do not.

Mixed effect model

I’m not going to run it becuase it takes a long time.

mixed.model = lmer(evaporation  ~ exchange + 
                     patch_size + 
                     disturbance  + 
                     exchange*patch_size + 
                     exchange*disturbance + 
                     patch_size*disturbance + 
                     exchange*patch_size*disturbance + 
                     (exchange||culture_ID) + 
                     (patch_size||culture_ID) + 
                     (disturbance||culture_ID) + 
                     (exchange*patch_size||culture_ID) + 
                     (exchange*disturbance||culture_ID) + 
                     (patch_size*disturbance||culture_ID) + 
                     (exchange*patch_size*disturbance||culture_ID), 
                   data=ds_for_evaporation,
                   REML = FALSE)

null.model = lm(evaporation ~
                  1, 
                data = ds_for_evaporation)

anova(mixed.model, null.model)

Other

Mixed effects models

To build the mixed effect models we will use the R package lme4. See page 6 of this PDF to know more about the syntaxis of this package.
To do model diagnostics of mixed effect models, I’m going to look at the following two plots (as suggested by Zuur et al. (2009), page 487):
- Quantile-quantile plots
- Partial residual plots
The effect size of the explaining variables is calculated in the mixed effect models as marginal and conditional r squared. The marginal r squared is how much variance is explained by the fixed effects. The conditional r squared is how much variance is explained by the fixed and the random effects. The marginal and conditional r squared are calculated using the package MuMIn. The computation is based on the methods of Nakagawa, Johnson, and Schielzeth (2017). For the coding and interpretation of these r squared check the documentation for the r.squaredGLMM function

Model selection

I am starting from the largest model and then simplifying because … (see statistical modelling course at ETH).
I am going to select the best model according to AIC. BIC is better for understanding and AIC for predicting. Halsey (2019) also suggests this approach instead of p values. I’m going to use AIC because I’m interested in knowign how much meta-ecosystem type contributed to the overall regional biomass.

Error bars

…

Body size

Aim

Data manipulation

Import

culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "morphology", "t0.RData"));t0 = morph_mvt
load(here("data", "morphology", "t1.RData"));t1 = morph_mvt
load(here("data", "morphology", "t2.RData"));t2 = morph_mvt
load(here("data", "morphology", "t3.RData"));t3 = morph_mvt
load(here("data", "morphology", "t4.RData"));t4 = morph_mvt
load(here("data", "morphology", "t5.RData"));t5 = morph_mvt
load(here("data", "morphology", "t6.RData"));t6 = morph_mvt
load(here("data", "morphology", "t7.RData"));t7 = morph_mvt
rm(morph_mvt)

Tidy

### --- Tidy t0 - t7 data-sets --- ###

#Column: time
t0$time = NA
t1$time = NA

#Column: replicate_video
t0$replicate_video[t0$file == "sample_00001"] = 1
t0$replicate_video[t0$file == "sample_00002"] = 2
t0$replicate_video[t0$file == "sample_00003"] = 3
t0$replicate_video[t0$file == "sample_00004"] = 4
t0$replicate_video[t0$file == "sample_00005"] = 5
t0$replicate_video[t0$file == "sample_00006"] = 6
t0$replicate_video[t0$file == "sample_00007"] = 7
t0$replicate_video[t0$file == "sample_00008"] = 8
t0$replicate_video[t0$file == "sample_00009"] = 9
t0$replicate_video[t0$file == "sample_00010"] = 10
t0$replicate_video[t0$file == "sample_00011"] = 11
t0$replicate_video[t0$file == "sample_00012"] = 12
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>% rename(replicate_video = replicate)
t7 = t7 %>% rename(replicate_video = replicate)


### --- Create ds dataset --- ###

long_t0 = t0 %>% slice(rep(1:n(), max(culture_info$culture_ID)))
ID_vector = NULL
ID_vector_elongating = NULL
for (ID in 1:max(culture_info$culture_ID)){
  ID_vector = rep(ID, times = nrow(t0))
  ID_vector_elongating = c(ID_vector_elongating, ID_vector)
}
long_t0$culture_ID = ID_vector_elongating
t0 = merge(culture_info,long_t0, by="culture_ID"); rm(long_t0)
t1 = merge(culture_info,t1,by="culture_ID")
t2 = merge(culture_info,t2,by="culture_ID")
t3 = merge(culture_info,t3,by="culture_ID")
t4 = merge(culture_info,t4,by="culture_ID")
t5 = merge(culture_info,t5,by="culture_ID")
t6 = merge(culture_info,t6,by="culture_ID")
t7 = merge(culture_info,t7,by="culture_ID")
ds = rbind(t0, t1, t2, t3, t4, t5, t6, t7); rm(t0, t1, t2, t3, t4, t5, t6, t7)

### --- Tidy ds data-set --- ###

#Column: day
ds$day = ds$time_point;
ds$day[ds$day=="t0"] = "0"
ds$day[ds$day=="t1"] = "4"
ds$day[ds$day=="t2"] = "8"
ds$day[ds$day=="t3"] = "12"
ds$day[ds$day=="t4"] = "16"
ds$day[ds$day=="t5"] = "20"
ds$day[ds$day=="t6"] = "24"
ds$day[ds$day=="t7"] = "28"
ds$day = as.numeric(ds$day)

#Column: time point
ds$time_point[ds$time_point=="t0"] = 0
ds$time_point[ds$time_point=="t1"] = 1
ds$time_point[ds$time_point=="t2"] = 2
ds$time_point[ds$time_point=="t3"] = 3
ds$time_point[ds$time_point=="t4"] = 4
ds$time_point[ds$time_point=="t5"] = 5
ds$time_point[ds$time_point=="t6"] = 6
ds$time_point[ds$time_point=="t7"] = 7
ds$time_point = as.character(ds$time_point)

#Column: eco_metaeco_type
ds$eco_metaeco_type = factor(ds$eco_metaeco_type, 
                             levels=c('S', 'S (S_S)', 'S (S_L)', 'M', 'M (M_M)', 'L', 'L (L_L)', 'L (S_L)'))

#Select useful columns
ds = ds %>% 
  select(culture_ID, 
         patch_size, 
         disturbance, 
         metaecosystem_type, 
         mean_area, 
         replicate_video, 
         day, 
         metaecosystem, 
         system_nr, 
         eco_metaeco_type)

#Reorder columns
ds = ds[, c("culture_ID", 
            "system_nr", 
            "disturbance", 
            "day",
            "patch_size", 
            "metaecosystem", 
            "metaecosystem_type", 
            "eco_metaeco_type", 
            "replicate_video",
            "mean_area")]

Create size classes data-set

As in Jacquet, Gounand, and Altermatt (2020) I will create 12 size classes.

#### --- PARAMETERS & INITIALISATION --- ###

nr_of_size_classes = 12
largest_size = max(ds$mean_area)
size_class_width = largest_size/nr_of_size_classes
size_class = NULL

### --- CREATE DATASET --- ###

size_class_boundaries = seq(0, largest_size, by = size_class_width)

for (class in 1:nr_of_size_classes){
  
  bin_lower_limit = size_class_boundaries[class]
  bin_upper_limit = size_class_boundaries[class+1]
  size_input = (size_class_boundaries[class] + size_class_boundaries[class + 1])/2
  
  size_class[[class]] = ds%>%
    filter(bin_lower_limit <= mean_area) %>%
    filter(mean_area <= bin_upper_limit) %>%
    group_by(culture_ID, 
             system_nr, 
             disturbance, 
             day, 
             patch_size, 
             metaecosystem, 
             metaecosystem_type, 
             eco_metaeco_type, 
             replicate_video) %>% #Group by video
    summarise(mean_abundance_across_videos = n()) %>%
    group_by(culture_ID, 
             system_nr, 
             disturbance, 
             day, 
             patch_size, 
             metaecosystem, 
             metaecosystem_type, 
             eco_metaeco_type) %>% #Group by ID
    summarise(abundance = mean(mean_abundance_across_videos)) %>%
    mutate(log_abundance = log(abundance)) %>%
    mutate(size_class = class) %>%
    mutate(size = size_input) %>%
    mutate(log_size = log(size))
  
}

ds_classes = rbind(size_class[[1]], size_class[[2]], size_class[[3]], size_class[[4]],
                  size_class[[5]], size_class[[6]], size_class[[7]], size_class[[8]],
                  size_class[[9]], size_class[[10]], size_class[[11]], size_class[[12]],)

Plot

Comparison of the three small patches

The code for plotting the size class distribution across small patches doesn’t work for some reason. It displays all plots as if they were the same.

Comparison of different time points

## [1] "Disturbance = low"

## [1] "Disturbance = high"

Single size class over time

Figures

Regional bioarea density (mean bioarea density between two patches) in meta-ecosystems of the same total area, but whose two patches have either the same size (red, medium-medium meta-ecosystem) or that have a smaller and larger patch (blue, small-large meta-ecosystem). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment.

Local biomass production in patches that are either connnected to a patch of the same size (green, small patches of small-small meta-ecosystemes) or to a patch of larger size (orange, small patches of small-large meta-ecosystems). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment. Problem: the small patches in green actually have different sample size than the one in orange. Does it still make sense to include the standard deviation?

Running time

## Time difference of 29.68097 secs

Bibliography

Halsey, Lewis G. 2019. “The reign of the p-value is over: What alternative analyses could we employ to fill the power vacuum?” Biology Letters 15 (5). https://doi.org/10.1098/rsbl.2019.0174.

Jacquet, Claire, Isabelle Gounand, and Florian Altermatt. 2020. “How pulse disturbances shape size-abundance pyramids.” Ecology Letters 23 (6): 1014–23. https://doi.org/10.1111/ele.13508.

Nakagawa, Shinichi, Paul C. D. Johnson, and Holger Schielzeth. 2017. “The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded.” Journal of the Royal Society Interface 14 (134). https://doi.org/10.1098/rsif.2017.0213.

Zuur, Alain F., Elena N. Ieno, Neil Walker, Anatoly A. Saveliev, and Graham M. Smith. 2009. Mixed effects models and extensions in ecology with R. Vol. 36. Statistics for Biology and Health. New York, NY: Springer New York. https://doi.org/10.1007/978-0-387-87458-6.

PatchSizePilot

Emanuele Giacomuzzo

2022-08-03

PatchSizePilot

Biomass

Aim

Data manipulation

Import

Tidy

Create regional data set

Evaporation

Tidy

Plot

Mixed effect model

Other

Mixed effects models

Model selection

Error bars

Body size

Aim

Data manipulation

Import

Tidy

Create size classes data-set

Plot

Comparison of the three small patches

Comparison of different time points

Single size class over time

Figures

Running time

Bibliography

PatchSizePilot

Emanuele Giacomuzzo

2022-08-03

PatchSizePilot

Biomass

Aim

Data manipulation

Import

Tidy

Create regional data set

Regional biomass

Medium-Medium vs Small-Large

Plot

Time as random effect

Tidy

Model selection

Time fitting

Time as fixed effect

Tidy

Model selection

Selected models

Local biomass

Small patches

Large patches

Evaporation

Tidy

Plot

Mixed effect model

Other

Mixed effects models

Model selection

Error bars

Body size

Aim

Data manipulation

Import

Tidy

Create size classes data-set

Plot

Comparison of the three small patches

Comparison of different time points

Single size class over time

Figures

Running time

Bibliography